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ABSTRACT 


A numerical method to solve two-phase turbulent flows with charged droplets 
in electrostatic field is presented. The ensemble-averaged Navier-Stokes equations 
and the electrostatic potential equation are solved using a finite volume method. 
The transitional turbulence field is described using multiple-time-scale turbulence 
equations. The equations of motion of droplets are solved using a Lagrangian 
particle tracking scheme, and the inter-phase momentum exchange is described by 
the Particle-In-Cell scheme. The electrostatic force caused by an applied electrical 
potential is calculated using the electrostatic field obtained by solving a Laplacian 
equation and the force exerted by charged droplets is calculated using the 
Coulombic force equation. The method is applied to solve electro-hydrodynamic 
sprays. The calculated droplet velocity distributions for droplet dispersions 
occurring in a stagnant surrounding are in good agreement with the measured data. 
For droplet dispersions occurring in a two-phase flow, the droplet trajectories are 
influenced by aerodynamic forces, the Coulombic force, and the applied 
electrostatic potential field. 
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NOMENCLATURE 


Ai p coefficient for i - th velocity component at a grid point p 

b vector of body force (= {b x ,b r }) 

Cp eddy viscosity coefficient 

Cpf constant coefficient (= 0 . 09 ) 

di diameter of a droplet 

Fi inter-phase force [Kg/(m 2 • sec 2 )], i = {jc, r } 
f a aerodynamic force 

f c Coulombic force 

f ES electrostatic force caused by applied electrical field 

h distance from meniscus to collecting surface 
k turbulent kinetic energy (= k p + k t ) 

k p turbulent kinetic energy in production range 

k t turbulent kinetic energy in dissipation range 

rri£ mass of a droplet 

N c i number of particles or droplets 
P r production rate 

p pressure 

p * initial guess for pressure 

p** corrected pressure 

p incremental pressure 

qj charge of a droplet 

R e £ droplet Reynolds number, R e £ - p\\ - \ pjdg / p 
r£ radius of a droplet 

external free stream velocity 
U 0 initial velocity of air jet 

uj initial guess for velocity 
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wy* predicted velocity 

uj** corrected velocity 

v ensemble averaged, gas phase velocity vector, v = {Uj\j=l,2} = {u, v} 

v £ velocity vector for a droplet, v £ = {u£ j |;=1,2} = {u # , V £ } 

x spatial coordinates (= {x, r}) 

At time-step size for numerical integration 
£(p permittivity of free space 

e p energy transfer rate 

£ t dissipation rate 

(f> electrostatic potential field 

ji molecular viscosity 

ji e effective viscosity (= fi + }X t ) 

H t turbulent viscosity 

p density of gas phase 


Superscripts 

n time-level for unsteady cases 

Subscripts 

e external boundary 

/ continuous gas phase 

t liquid droplets 

o initial conditions 

O average value 
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1. Introduction 

In electrohydrodynamic (EHD) sprays, dispersion of particles or droplets is 
caused by electrical force and aerodynamic force. Application of EHD sprays can 
be found in spray paintings, aerosols, fabrication of solar panels, and dust 
collecting mechanisms. The use of electrical force also offers the possibility of 
controlling dispersion of liquid fuel droplets in combustors and aerospace 
propulsion systems. A number of experimental and theoretical investigations have 
been made during past decades on two-phase flows for applications in combustors 
and aerospace propulsion systems. However, only a limited amount of research 
efforts have been made for two-phase flows with a charged discrete phase. 
Certainly, more experimental and theoretical research efforts will be made in the 
future to optimize dispersion of particles or droplets in various engineering 
applications. A few experimental and theoretical works closely related to the 
present work are discussed below. 

Experimental data on electrostatic spray emitted from a meniscus used in the 
present work can be found in ref. [1]. Subsequently, the experimental case was 
studied numerically by Ganan-Calvo et al. [2] using a Lagrangian particle tracking 
technique. The charged droplets emanating from the meniscus are accelerated by 
an applied electric potential and by an air jet surrounding the meniscus [3]. 
However, detailed data on the air jet are not available. In the work of Ganan-Calvo 
et al., the aerodynamic force acting on a droplet was calculated assuming that the 
surrounding air is in a stagnant state. Numerical droplets are seeded at locations 
60-80 droplet radii downstream of the conical meniscus. Their calculated velocity 
profiles are in good agreement with measured data. The good comparison between 
the calculated results and the measured droplet velocity profiles indicates that the 
influence of the airflow on the droplet motion is negligibly small. Experimental 
investigations on dispersion of chemically reacting, charged droplets have also 
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appeared in recent years [4]. However, detailed measured data for two-phase 
turbulent flows with charged discrete phase are still very scarce. 

In gas-droplet two-phase flows, inter-phase mass, momentum, and energy 
exchanges exist. For the two-phase flow considered herein, the influence of 
evaporation of liquid droplets on the gas phase is negligible since the time interval 
that droplets are exposed to airflow is very small. Also, energy exchange is ignored 
since the fluid flow is in an isothermal state. In the present numerical investigation, 
the inter-phase momentum exchange is resolved using the particle-in-cell (PIC) 
scheme of Crowe et al. [5]. Implementation of the PIC is described in more detail 
in the "Numerical methods" section. 

The numerical method used in this work is a pressure-based Navier-Stokes 
equations solver that incorporates a pressure-staggered mesh and an incremental 
pressure equation for the conservation of mass. The accuracy of the numerical 
method has been validated by solving a number of flow cases. Calculations of two- 
dimensional self-sustained oscillatory flows over a circular cylinder and a square 
cylinder can be found in ref. [6], and calculations of three-dimensional lid-driven 
cavity flow and a laminar flow through a curved duct can be found in ref. [7]. 
Further application of the numerical method for various laminar flows, complex 
turbulent flows, incompressible flows, compressible flows, steady flows, unsteady 
flows, and chemically reacting flows can be found in refs. [8-11] and the references 
cited therein. 

It has been shown that the multiple-time-scale (M-S) turbulence equations used 
in this work yield accurate numerical results for various complex turbulent flows. 
The numerical results for various incompressible and compressible turbulent flows 
(e.g., an unsteady transitional flow over an oscillating airfoil, circular jets in 
crossflow, transonic flows with shock wave - boundary layer interactions, a 
confined swirling jet, and divergent channel flows) obtained using the M-S 
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equations are in as good agreement with the measured data as those obtained using 
optimized k-8 turbulence models, algebraic Reynolds stress turbulence models 
(ARSM), or Reynolds stress turbulence models (RSM) for each flow case. It has 
also been shown that the M-S equations can resolve ignition delay, flame 
thickness, and chemical reaction - turbulence interaction occurring in a 
compressible shear layer [8]. It can be found in ref. [8] that the distribution of 
chemical species obtained using the M-S equations are in much closer agreement 
with the measured data than those obtained using k-e turbulence equations or a 
Monte-Carlo probability density function method. The physical aspects of the M-S 
equations are briefly discussed in this paper. 

2. Theoretical Analysis 

The mathematical equations governing an unsteady, two-phase, turbulent flows 
with a charged discrete phase include: the Navier-Stokes equations, turbulence 
equations, Lagrangian equations of motion for discrete phase, the Coulombic force 
equation, and an electrostatic potential equation. These equations are discussed 
below. 

Conservation equations for gas phase flow 

The electrostatic force acting on a droplet is contributed by the applied 
electrostatic potential and by charged droplets. The spatial location of each charged 
droplet is a function of time and, hence, a time-accurate solution technique needs to 
be used to solve two-phase turbulent flows with charged droplets. The ensemble- 
averaged conservation of mass equation is given as 

+ ; |M) = 0. (1) 

where p is the density and {w, v} are the ensemble-averaged velocities. 
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The ensemble averaged conservation of linear momentum equations in 
axisymmetric coordinates are given as 
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where p is the ensemble-averaged pressure, {F x , F r } is the inter-phase force, 
p e = p + p t is the effective viscosity, and p and p t are the molecular and 
turbulent viscosities, respectively. A numerical procedure for the inter-phase 
momentum exchange term is described later in the "Numerical methods" section. 


Multiple-time-scale turbulence equations 

The M-S equations are based on a simplified split-spectrum method [10]. In the 

method, the turbulent kinetic energy is split into turbulent kinetic energy in the low 
frequency range (k p ) and that in the high frequency range (k t ). Recall that the 

energy-containing eddies are generated by the instability of the mean fluid flow, 
the energy-containing eddies cascade to finer eddies, and the fine scale eddies are 

dissipated by the viscous force. Hence, the energy containing eddies are 
characterized by low frequencies and large values of k p /k t and the fine scale 

eddies are characterized by higher frequencies and small values of k p /k t . The 

capability to resolve the cascade is achieved by solving the convection-diffusion 

equations for the spectrum-split turbulent kinetic energies. Obviously, single-time- 
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scale turbulence models can not resolve the cascade of turbulent kinetic energy 
since spectrum- split turbulence quantities are not solved for in these turbulence 
models [10]. 

In complex turbulent flows, the production rate (P r ) and dissipation rate (e t ) of 
the turbulent kinetic energy vary widely in space so that the shape and the 
frequency domain of the spectral density also vary widely in space. Such a state of 
turbulence is called "nonequilibrium turbulence." The capability to resolve the 

nonequilibrium turbulence phenomena originates from describing the turbulence 
length scale and the turbulent viscosity using the energy transfer rate ( e p ). In the 

M-S equations [10], the eddy viscosity coefficient that depends on the strength of 
nonequilibrium turbulence is given as; 

Cpi = Cjjf £ t / e p (4) 

where is a constant coefficient. A few experimental and theoretical 
investigations on the dependence of the eddy viscosity coefficient on P r /e t can be 
found in refs. [12-14]. In single-time-scale turbulence models, the turbulence 
length scale is defined using the dissipation rate and, hence, these turbulence 
models can not resolve the nonequilibrium turbulence phenomena. 


Lagrangian description of droplet equations of motion 

The Lagrangian equations of motion for discrete phase are given as [5, 15-17]; 


dt 


= V£ 





(5) 

( 6 ) 


9 



where trig is the mass of a droplet, =— c d p\v- v ^|(v- v^);rr/ +b is the sum of 

the aerodynamic force and the electrostatic force, nr} is the wetted area of a 
droplet, and b is the electrostatic force. The drag coefficient is given as 


c d = 


24 

^ el 




(7) 


where the droplet Reynolds number is defined as 
_p|v-v/ 


R. 
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The body force acting on a charged droplet in electrostatic potential field is 
given as 
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where qj is the charge of a droplet, (f> is the electrostatic potential field, Nj is the 
number of charged particles or droplets in the flow field, is the permittivity of 
free space, R j k = xj - x k and Rj k = R j k . For the class of problems considered in 

this study, the applied electrostatic field is in a steady state and the governing 
partial differential equation is given as [18]; 


V 2 </S = 0. 


( 10 ) 


3 Numerical Method 


The unsteady gas phase equations are solved using a finite volume method. In 
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the method, the velocities and turbulence quantities are located at the grid points 
and the pressure and the electrostatic potential are located at the centroid of a cell 
formed by the four adjacent velocity grid points. Interaction between the 
isothermal continuous gas phase and the discrete droplet phase are accounted for 
by considering the inter-phase exchanges of momentum. The transient solution for 
the gas phase is obtained by solving the flow equations iteratively at each time- 
level. The numerical method for the gas phase equations that include the inter- 
phase source term, particle equations of motion, and the applied electrostatic 
potential field are described below. 


Numerical method for gas phase equations 

The numerical method is a pressure based Navier-Stokes equations solver in 
which the predicted velocity is obtained by solving the momentum equations and 
the divergence free, corrected velocity is obtained by solving an incremental 
pressure equation derived from the conservation of mass equation. 

Let u* and p* be the initial guesses for the velocity and pressure for a new 

time-level, respectively. Applying a finite volume method [6-11] to the momentum 
equations yields 


(pC\ + \ p - X Aij&ijk ~ ~T~ A V + pC 2 u”p + S t + Spi (1 1 ) 
&=1 


where Cj = C 2 = f... -p—dx, AV is the volume of a cell, At * = t n ~t n 1 is the 
Ay A tf J 

time-step size, n denotes the time-level, nb is the number of neighboring grid 
points, A i p is determined from the power-law upwind differencing scheme, 5/ 

represents source terms contributed by non-orthogonal mesh and other velocity 

* 

components, and represents the inter-phase source term. The predicted 
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velocity, w**, is obtained by solving eq. (11). 

The inter-phase momentum transfer from the discrete phase toward the gas 
phase is contributed by gain or loss of momentum of the discrete phase. In the 
context of the PIC scheme, the inter-phase source term for the discrete momentum 
equation for each cell is obtained by summing up the source terms contributed by 
all the droplets inside the cell. 



=-Lf'-rj F t dxdt 

A tfJt=t n J AV 1 




( 12 ) 


where iqj is the droplet velocity in the i - th coordinate direction. 

The predicted velocity field does not satisfy the conservation of mass until the 
solutions are fully converged. The corrected velocity and pressure that satisfy the 
conservation of mass and the momentum equations can be written as 


«/ = 


Mf I 




( 13 ) 


Then the corrected discrete momentum equation can be written as 
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The inter-phase source term is not updated during the iterative solution process of 
the gas phase equations. Therefore, the time-step size of the gas phase for two- 
phase unsteady flows needs to be much smaller than that for single phase unsteady 
flows. Otherwise, the inter-phase momentum transfer may not be resolved 
accurately. Subtracting eq. (11) from eq. (14) yields the relationship between the 
incremental velocity and the incremental pressure given as 


u 


/ 

i 



(15) 


where A u . = A v/(pCi +A* p ), and the contributions made by neighboring grid 

points have been disregarded. However, this simplification does not incur any mass 
imbalance in the converged solution since the incremental pressure is driven only 
by the mass imbalance as shown below. 

The incremental pressure equation is obtained by inserting eq. (13) into eq. (1) 
and is given as, after some rearrangement. 



where the right hand side represents the mass imbalance and is identical to the 
original conservation of mass equation. Therefore, the incremental pressure is 
driven only by the mass imbalance and the simplifications introduced during 
derivation of the incremental pressure equation do not incur any mass imbalance in 
the converged solution. 

Applying the finite volume method to eq. (16) yields 

ApP p ILAkPk y P Mj fijds (17) 
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where ds is the boundary surface of a pressure control volume. The incremental 
pressure is obtained by solving eq. (17), the incremental velocity is calculated 
using eq. (15), and the corrected velocity and pressure are calculated using eq. (13). 

The turbulence equations are solved by the same finite volume method. At each 
time-level, the conservation of momentum equation, the incremental pressure 
equation, and the turbulence equations are solved iteratively until the velocity and 
pressure satisfy the prescribed convergence criterion. In principle, the incremental 
pressure equation needs to be solved iteratively until the residual term vanishes 
and, in such a case, the converged solution satisfies the conservation of mass and 
momentum equations exactly within the context of the difference approximation. 


Particle Equations of Motion 

The particle equations of motion are integrated using the Crank-Nicholson 
scheme. Let A t# be the time-step size for the discrete phase. Then eqs. (5) and (6) 

can be written as; 


*?-*? 1 _ 1 v n+l/2 

A t t ~2 l 


mi 


v ?“ v ? 1 _l F n+l/2 


( 18 ) 

(19) 


Eqs. (18) and (19) form a set of nonlinear equations since depends on the 
droplet velocity as well as the gas-phase velocity that is a function of the spatial 
location. Ideally, eqs. (18) and (19) needs to be solved iteratively. Instead, a very 
small time-step size can be used to semi-implicitly integrate eqs. (18) and (19). The 
particle velocity and the force on the right hand side of eqs. (18) and (19) are 
calculated using 
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( 20 ) 


v? +1/2 4( V ?-‘+v?) 

F " +1/2 = I ( F " _1 + F ") 

where \g and F" are evaluated using an explicit time-integration scheme. In gas- 
liquid two-phase flows, the physical time-scales as well as the numerically stable 
time-step-sizes for the continuous and the discrete phases are quite different. The 
use of a very small time-step size is necessary to accurately resolve the inter-phase 
momentum transfer. A very small time-step size also needs to be used for 
numerical stability. The use of a time-step size approximately two orders of 
magnitude smaller than that of the continuous phase yields stable and strongly 
convergent results. A similar conclusion can also be found in Raju and Sirignano 
[17]. 

The electrostatic force acting on a charged droplet is caused by the applied 
electrostatic potential field and by the Coulombic force acting among charged 
droplets. The force exerted by the applied electrical field is calculated using the 
electrical potential field obtained by solving a Laplacian equation. In numerical 
calculation, the electrical potential is defined at the center of a cell and the 
Laplacian equation is solved using the same procedure as that for the incremental 
pressure equation [6-1 1], The force exerted by charged droplets is calculated 
directly using the Coulombic force equation. 

The liquid spray is represented by a finite number of droplets with different size 
groups. Let Mg ^ be the mass flow rate of the droplets in the k - th size group. 

Then the total mass flow rate. Mg, of the droplets is given as 



^siz . 

£ M £,k 

*=i 


( 21 ) 
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where N s i z represents the number of all droplet size groups. The number flow rate 
of the physical droplets in the k — th size group is given as 


*lk ~ 4 3 


( 22 ) 


The number of physical droplets for each size group to be introduced at each time- 
level of the discrete phase is given as 

N$ = rj k Atg + Resid (wf 1 ) (23) 

where Atg is the time-step size for the discrete phase, and Resid|^ -1 j is the 

residual carried over from the previous time-level of droplet injection. 

The charge carried by each droplet is calculated using 


qj=q 


Ki|“ 

\dg j 


(24) 


where dg is the average diameter of the droplets. The charge carried by an average 
size droplet, q, is obtained by dividing the measured total current by the number 

flow rate of average-sized droplets. According to the experimental work of Gomez 
and Tang [19], a ~ 2 for dg smaller them approximately 50 }Jm, and a ~ 1.5 for 

substantially larger droplets. Theoretical investigations on charge-to-mass ratio 
show that a ~ 1.5 [20]. In the present work, a ~ 1.5 is used following the work of 
Ganan-Calvo et al. [2]. 


4. Results 


Two different cases of droplet dispersion are considered herein. In the first case, 
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calculations are made for dispersions of charged droplets occurring in a stagnant 
surrounding and the results are compared with available numerical results and 
measured data [1,2]. In the second case, dispersion of charged droplets coupled 
with fluid flow is calculated and the capability of the numerical method to solve 
droplet or particle dispersions occurring in laminar and turbulent flows is 
demonstrated. The calculated results are analyzed to investigate the trajectory of a 
droplet and the forces acting on the droplet. Physical data used in both calculations 
are described below. 

The electrostatic spray emanating from a meniscus and the extent of the 
computational domain are shown schematically in Fig. 1. The length of the needle 
is 0.028 m, the inner radius (h) is 0.00025 m, and the outer radius is 0.00050 m. 
The distance between the meniscus and the collecting surface ( H ) is 0.0254 m. 

The other physical dimensions of the computational domain are: a = 0.014 m, 

b = 1.24 H, c = H, and d = 0.95 H. The applied electric potential difference is 3040 
volts, and the current carried by charged droplets is 4.3 x 10 -8 Ampere. The liquid 

is heptane doped with an antistatic additive (ST ADIS 450, Du Pont). The liquid 
density is 685 Kg/m 3 , and the volumetric flow rate is 2.3 x 10 -9 m 3 / sec. The 
measured data show that droplet diameters are distributed in between 35 ~ 50 pm 

[2]. In the present work, the spray is represented by 5 different size groups with the 
diameters distributed between dg - g s d# and dg + cr s d£, where dp =40 pm is the 
mean diameter and <J S = 0.20 is the standard deviation. The number flow rate of 

each size droplet is calculated using the Gaussian distribution. In this case, the 
charge carried by an average-sized droplet is 6.26 x 10 -13 Coul. The velocity of the 

liquid spray leaving the needle is 0.012 ml sec. The droplets are accelerated by the 
applied electric potential and by an air jet surrounding the needle [3]. However, 
detailed data for the air jet are not available. The annulus of the air jet is assumed 
to be 0.00025 m as shown in Fig. 1. 
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Dispersion of droplets in stagnant surrounding 

The electrostatic spray emanating from a conical meniscus considered herein 
can be found in refs. [1,2]. Following the previous work [2], it is assumed that the 
airflow is negligibly small. The mesh for the entire domain and that near the 
injector are shown in Fig. 2.(a) and Fig. 2.(b), respectively. The smallest mesh size 
is 0.5 x 10 -4 m and the mesh away from the injector is stretched by a factor of 1.25. 
Note that the radius of the liquid injector is 0.25 x 10 _3 m while that for the 
collector surface is 0.254 x 10 _1 m. The length scale of the injector is two orders of 
magnitude smaller than that of the entire domain as well as the droplet collector 
surface. The use of such a large stretch ratio is necessary to discretize a 
computational domain characterized by largely disparate length scales. 

The boundary conditions for the electrical potential field are: 0 = 3040 volts for 
the meniscus, 0 = 0 along the surface of the droplet collector, and d<p/dn = 0 along 
the center line and the entire external boundary. The air jet and air duct are treated 
as a free space so that the calculated potential field can be compared directly with 
the analytical and numerical results presented in ref. [2]. 

The calculated electrical potential field is shown in Fig. 3(a). The potential 
along the centerline obtained in the present study is in good comparison with 
analytical and numerical results [2] as shown in Fig. 3(b). The good comparison 
indicates that the present numerical method can resolve the highly steep and almost 
singular potential field near the meniscus. The calculated lines of force and the 
magnitude of force are shown in Figs. 3(c) and 3(d), respectively. Fig. 3(d) also 
shows that a strong electrical force field is concentrated in the region very close to 
the conical meniscus. 

The number flow rate of average-sized droplets is calculated by dividing the 
volumetric flow rate of the conducting liquid by the volume of an average-sized 
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droplet, and the charge carried by an average-sized droplet is obtained by dividing 
the current by the number flow rate. Thus the number flow rate is inversely 
proportional to the cubic power of the droplet diameter, while the charge carried by 
an average-sized droplet is proportional to the cubic power of the diameter. The 
lines of force are closely aligned with the centerline in the spray region. Hence, the 
lateral dispersion of droplets is mostly caused by the Coulombic force acting 
among charged droplets and, hence, the calculated droplet dispersion in the radial 
direction depends strongly on the average droplet size. A few trial calculations 
revealed that the droplet dispersion also depends strongly on the mode of 
numerical droplet injection. The cause for such a dependence is discussed in the 
following "Two-phase flow with charged droplets" section. For the reasons 
discussed above, small uncertainties in measured current, volumetric flow rate, and 
the average droplet size can produce significantly different numerical results. To 
reduce uncertainty, the numerical droplets are seeded at x - 0.012 m in the present 

calculation. The lateral seed location is obtained using a random number generator, 
and the initial droplet velocity is obtained from measured data at the seed location. 

The time-step size for numerical integration of droplet equations of motion is 
A t£ = At y / 100, where A tf- 0.175 x 10" 4 sec is the time-step size for the gas 

phase to be discussed in the following section. Calculations were also made using 
At £ = At f / 200 to confirm the independence of the calculated results on the time- 

step size. The numerical results obtained using two different time-step sizes are 
practically identical. The initial droplet velocities and the calculated results at 
downstream locations are shown in Fig. 4. It can be found in the figure that the 
droplet velocity profiles obtained in the present study are in good agreement with 
the measured data. However, the present results exhibit wider later dispersion of 
droplets than those of ref. [2], and a few droplets disperse beyond r = 0.01 m. 
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Two-phase flow with charged droplets 

The inlet boundary of the air jet is located at the tip of the meniscus. The initial 
velocity of the air jet is 22 m/sec along the 45-degree inclined meniscus surface. 
The mesh and the applied electrical potential are the same as those described in the 

"Dispersion of droplets in stagnant surrounding" section. The time-step size for the 
fluid flow is Atf=T c /N t , where T c = 0.0035 sec is a characteristic time for a 

droplet to reach the surface of the collector and N t = 200 is the number of time- 

steps per characteristic time. The boundary conditions are as follows. At the 
upstream boundary above the air duct, u-U e , v = 0, and a weak turbulence field 

is prescribed. Along the inclined surface of the liquid jet, the axial velocity of air 
is set equal to the liquid injection velocity, and v = k p =k t =0 and 

d£pjdn = de t /dn = 0. Along the center line and at the exit boundary, v = 0 and a 

vanishing gradient boundary condition is used for all other variables. At the 
external far field boundary, u-JJ e , v = 0, and a vanishing gradient boundary 

condition is used for all other variables. Along the solid wall boundary of the 
droplet collecting surface, u = v = k p =k t = 0 and depjdn = de t /dn = 0. External 

far fields are in a stagnant state in many EHD sprays. For the two-phase flow case 
considered in this section, more than 90 per cent of the domain is in a stagnant 
state. A stagnant state is certainly a trivial solution of the Navier-Stokes equations. 
However, numerically resolving such a stagnant state that extends through a large 
portion of the flow domain is not a trivial matter. To overcome numerical 
instability caused by the no fluid motion in a large extent of the domain, a small 
amount of free stream velocity, U e = 0.035 U 0 , has been used. Since the fluid flow 

and the particle dispersion occur in a small region along the center line, the 
external velocity does not significantly alter the two-phase fluid flow. Developing 
a numerical method to solve fluid flow occurring in a very small region of a large 
stagnant surrounding is a formidable task and it will be reported separately. 
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The calculated velocity vectors, pressure contours, and the turbulent kinetic 
energy contours for the gas-phase flow are shown in Figs. 5(a)~5(c), respectively. 
These figures show that a significant fluid motion occurs only in a narrow region 
along the center line. A very weak turbulence field is developed only along the air 
jet and the entire domain, including the center region, remains in a practically 
laminar state. 

The convergence history for velocities, pressure, and conservation of mass is 
shown in Fig. 6(a) and that for turbulence variables is shown in Fig. 6(b). The L2 
error norm for each flow variable is defined as 


lk(«)ll= 


n Y ~ 1 »v 1 


il/2 


X i l{{<T- a ij)/ a } /k- 2)0*3, -2) 


(25) 


L *=2 J = 2 


where n x and n y are the number of grid points in the axial and radial directions, 
respectively, and a = max = 2 ,n x - 1 and j = 2,n y - 1 j is a normalizing 

factor. The global conservation of mass error is defined as 

I e CM | = (***i n ~ ™out )A**m ( 26 ) 


where m- m and rh out represent mass fluxes entering and leaving the domain, 
respectively. The small external velocity does not completely suppress the velocity 
and pressure wiggles. Yet, the present numerical method yields strongly 
convergent results as shown in these figures. It is shown in Fig. 6(b) that the 
multiple-time-scale turbulence equations yield strongly converged results for the 
weak turbulence field developing in a flow field with largely disparate length 
scales. The capability of the turbulence equations to accurately resolve complex 
transitional turbulence field can be found in ref. [1 1]. 
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The droplets are seeded into the flow field after the flow field has reached a 
converged state. The droplet seed location confined within 
{(x, y) 1 0.0035 < x < 0.0042, 0. < y < 0.0007} is calculated using a random number 
generator and the initial velocity of droplets is 5 ml sec. As in the previous section, 
the use of At# = At f /100 yields numerical results that are independent of the time- 

step size. The calculated streaklines are shown in Fig. 7. The convergence history 
for flow and turbulence variables are shown in Figs. 8(a) and 8(b), respectively. 
The number of droplets residing in the flow field versus time-level is shown in Fig. 
8(c). As a steady state is approached, approximately 290 droplets exist in the flow 
field. For the two-phase flow case, the inter-phase momentum source term 
contributed by the discrete phase continuously disturbs the velocity field of the gas 
phase. To prevent accumulation of errors, the flow equations are solved iteratively 
for at least a prescribed number of iterations at each time-level. As shown in Figs 
8(a) and 8(b), the use of 5 iterations per time-level let the flow and turbulence 
variables stay within a consistently converged state. After the prescribed minimum 
number of iterations is performed, calculations advance to a new time-level if 
either a prescribed convergence criterion is met or the number of iterations exceeds 

the prescribed maximum number of iterations. The convergence criterion used is 
| e CM | < 1.0 x 10 -3 and the maximum number of iterations is 25 per time-level. 

The trajectory and velocity components of a droplet are shown in Fig. 9(a), and 
the axial and radial forces acting on the droplet are shown in Figs. 9(b) and 9(c), 
respectively. Immediately after the droplet is seeded, it is accelerated by the 
applied electrostatic force and aerodynamic force. Once the droplet attains its peak 
velocity, the aerodynamic force begins to decelerate the droplet. It can be found in 
Fig. 9(c) that the lateral dispersion of droplets is mostly caused by the Coulombic 
force. When a droplet is seeded at a relatively remote location from other droplets, 
the Coulombic force acting on the droplet is small. However, as soon as a new 
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droplet is seeded at a close location, the Coulombic force acting on the droplet is 
increased instantly. The dependence of droplet dispersion on the mode of 
numerical droplet seeding is caused by this phenomenon. As the particle travels 
toward the downstream direction, the Coulombic force becomes negligibly small 
due to the dispersion of droplets. At such downstream locations, a droplet can 
attain an equilibrium velocity that is balanced by the applied electrostatic force and 
the drag force. 

5. Conclusions and discussion 

A time-accurate numerical method to solve two-phase turbulent flows with 
charged particles or droplets in electrostatic field is presented. The method can 
solve dispersions of a discrete phase occurring in a stagnant surrounding, in 
laminar flows, and in turbulent flows. Therefore, the method can be used to 
investigate the dispersion of a charged discrete phase encountered in various 
engineering applications. The capability to analyze atomization of electrically 
conducting liquids is not included in this study. Atomization is an important , on- 
going, research subject by itself even for liquids that do not involve electrical force. 
The present method can be used advantageously for numerical investigation of 
such problems as spray painting and dust collecting mechanisms that do not 
involve atomization of electrically conducting liquids. For problems that involve 
atomization, such as aerosols, fabrication of solar panels, and the problem 
considered in the present study, the method provides limited capability to study 
dispersion of droplets that involve electrical force. 

For the droplet dispersion occurring in a stagnant surrounding, the calculated 
droplet velocity profiles are in good agreement with the measured data. The 
calculated results obtained in the present study exhibit wider lateral dispersion than 
that obtained by Ganan-Calvo et al. [2]. The lateral dispersion is mostly caused by 
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the Coulombic force acting among charged droplets. It is shown that the 
Coulombic force depends strongly on the charge carried by each droplet and the 
mode of numerical droplet injection. In ref. [2] and here, the charge carried by an 
average-sized droplet is calculated by dividing the current by the number flow rate 
of droplets. The number flow rate is inversely proportional to the cubic power of 
the average droplet diameter. Hence a small uncertainty in the measured droplet 
size, current, and volumetric flow rate can cause a substantial difference in 
calculated results. 

For the flow case considered in this study, a large extent of the domain is in a 
nearly stagnant state. It is found that the large stagnant surrounding causes 
numerical instability. Introducing a small external velocity suppresses the 
numerical instability and yields strongly convergent and physically correct results. 
A weak turbulence field is developed only near the exit of the air jet and the entire 
flow field remains in an almost laminar state. It is shown that the multiple-time- 
scale turbulence equations do not cause difficulty in obtaining a highly convergent 
result. In fact, the turbulence equations yield strongly converged results for the 
weak turbulence field developing in a flow field with largely disparate length 
scales. Comparisons of the aerodynamic force, electrostatic force caused by the 
applied electrical field, and the Coulombic force acting on a charged droplet help 
to better understand the characteristics of dispersion of charged particles or 
droplets occurring in two-phase flows. Measured data for two-phase flows with a 
charged discrete phase are very scarce, and more detailed measured data are 
certainly necessary for further improvement and verification of numerical methods 
to solve two-phase flows with a charged discrete phase. 

A few difficulties in numerical simulations of two-phase flows with a charged 
discrete phase are discussed below. Usually, the time-step size of the discrete phase 
in two-phase flows needs to be by far smaller than that for the gas phase in order to 
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obtain numerical results that are independent of the time-step size used. Otherwise, 
the inter-phase momentum transfer may not be resolved correctly. The time-step 
size of the discrete phase for two-phase flows with charged droplets needs to be 
even smaller than that without Coulombic force. If the time-step size for the 
discrete phase is not sufficiently small, then droplets can approach infinitely close 
to each other and these droplets may be subjected to a large Coulomb force since 
the Coulomb force between charged droplets is inversely proportional to the square 
of the distance between the droplets. In such a case, a large Coulomb force 
accelerates droplets so greatly that an unphysically large inter-phase momentum 
transfer rate can be produced and the numerical method for the gas-phase may fail 
to yield a converged solution. Only the use of a very small time-step size can yield 
a physically correct droplet distribution pattern. 
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1 . Nomenclature and computational domain for EHD spray. 
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(a) mesh for entire domain, 



(b) mesh near injector. 


2. 109x87 Mesh, 
29 



(a) contour plot, 0 < 0 < 3040 volts, and A0 = 160 



3. Applied electrostatic potential field 
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(c) Lines of force. 








(b) pressure contours, -10 <p<20 Newton/ m , A p = 1.50, 
5. Calculated flow field 
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2 2 

(c) turbulent kinetic energy contours, 0<fc<8m / sec , Ak = 0.8. 


Figure 5 - continued. 
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6. Convergence history for fluid flow, 
(a) flow variables, (b) turbulence variables. 
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8. Convergence history for two-phase flow, 

(a) flow variables, (b) turbulence variables, (c) number of droplets. 
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